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Abstract 

We present results from an extensive set of simulations of gravitational fragmentation in the 
presence of magnetic fields and ambipolar diffusion. The thin-sheet approximation is employed, 
with an ambient magnetic field that is oriented perpendicular to the plane of the sheet. Nonlinear 
development of fragmentation instability leads to substantial irregular structure and distribu- 
tions of fragment spacings, fragment masses, shapes, and velocity patterns in model clouds. We 
study the effect of dimensionless free parameters that characterize the initial mass-to-flux ratio, 
neutral-ion coupling, and external pressure associated with the sheet. The average fragmenta- 
tion spacing in the nonlinear phase of evolution is in excellent agreement with the prediction 
of linear perturbation theory. Both significantly subcritical and highly supercritical clouds have 
average fragmentation scales (A) ~ 2ttZo, where Zq is the initial half-thickness of the sheet. 
In contrast, the qualitatively unique transcritical modes can have (A) that is at least several 
times larger. Conversely, fragmentation dominated by external pressure can yield dense cluster 
formation with much smaller values of (A). The time scale for nonlinear growth and runaway 
collapse of the first core is ~ 10 times the calculated growth time T g , m of the eigenmode with 
minimum growth time, when starting from a uniform background state with small-amplitude 
white-noise perturbations. Subcritical and transcritical models typically evolve on a significantly 
longer time scale than the supercritical models. Infall motions in the nonlinear fully-developed 
contracting cores are subsonic on the core scale in subcritical and transcritical clouds, but are 
somewhat supersonic in supercritical clouds. Core mass distributions are sharply peaked with 
a steep decline to large masses, consistent with the existence of a preferred mass scale for each 
unique set of dimensionless free parameters. However, a sum total of results for various initial 
mass-to-flux ratios yields a broad distribution reminiscent of observed core mass distributions. 
Core shapes are mostly near-circular in the plane of the sheet for subcritical clouds, but be- 
come progressively more elongated for clouds with increasing initial mass-to-flux ratio. Field 
lines above the cloud midplane remain closest to vertical in the ambipolar-drift driven core 
formation in subcritical clouds, and there is increasing amount of magnetic field curvature for 
clouds of increasing mass-to-flux ratio. Based on our results, we conclude that fragmentation 
spacings, magnitude of infall motions, core shapes, and, especially, the curvature of magnetic 
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field morphology, may serve as indirect observational means of determining a cloud's ambient 
mass-to-flux ratio. 
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1. Introduction 



1.1. Molecular Cloud Cores 



Star formation occurs in dense cores within interstellar molecular clouds. The exis- 
tence and pro perties of dense cores are well established by studies of molecular sp ectral 
line emission ( Mvers fc Benson . 1983 : Benson fc Mversl . 1989 : jTiina et all 19991). sub- 
millim eter dust emission (jWard -Thompson et al. . 1994 : Andre ct al 



Kirk et al 



2005 ) , and infrared absorption (jBacmann et all 2000t Teixeira et al I l2005HLada et al 



20071 ). The core formation process has been the subject of intense theoretical study for the 



past few decades. Ideas range fro m uninhibited gravitational fragmentation instability 
( Jeansl 1929 : Larsonl 19851 20031). to ambipolar di f fusion in a magn etically supported 
cloud (|Mestel fc Spitzerl . Il956l iMouschoviasl ll978HShu et all Il987t), t o a very rapid 
fragmentation due t o pre-existing turbulent flows (jPadoan et al 1 11997k [Klessenl . l200l[ 



iGammie et al. , 2003 ) . The last two ideas are related to scenarios for the nature of rela- 
tively low-den sity molecular cloud env elopes, which contain most of the mass of molecu- 
lar clouds (see Goldsmith et al. . 2008). In one scenario, they are magnetically-dominated 
and evolve due to ambipolar diffusion, a gravita tionally-driven redistribution of mass 
amongst magnetic flux tubes ( Mouschovias . Il978t ) . The r elatively low g lobal star forma- 
tion rate and efficiency in Galactic molecular clouds (see IMcKee . 19891 ) is us ed to argue 



Shu et al 



that molecular cloud e nvelopes have a subcritical mass-to- flux ratio (see e.g 
1999; El megree nl l2007l) . Indirect empirical evidence based on the Chandrasekhar-Fermi 
(CF) method and velocity anisotr opy measurements does imply that the low density (n ~ 



100 cm 3 ) regions are subcritical (jCortes et all 120051 : iHever et all 120081 ). The contrasting 



view is that the mass-to-flux ratio is supercritical on large scales, and that the star for- 
matio n rate and efficiency are governed by supersonic turbulence ([Mac Low fc Klessen . 
20041) . This requires either a continual replenishment of rapidly dissipating turbulence, 



or a relatively rapid dispersal of the cloud envelopes. The focus of this paper is the frag- 
mentation of embedded and relatively quiescent dense regions, and not on the nature 
of the larger envelopes. We study the effect of small-amplitude perturbations on a vari- 
ety of cloud models with different mass-to-flux ratios, degrees of magnetic coupling, and 
external pressures. 

For the purpose of this paper, we refer to any individual unit of star formation, which 
leads typically to a single or small multiple star system, as simply a "core" . However, as- 
tronomers often subdivide this concept into two observational categories, that of "prestel- 
lar core" and "prestellar condensation" . The former term is often used to describe the 
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extended (mean density n ~ 10 5 cm -3 , size s ~ 0.1 pc, spacing A ~ 0.25 pc) objects 
in regions of distributed star formation like the Taurus molecular cloud, while the latter 
term typically describes the more compact (n > 10 6 — 10 7 cm -3 , s ~ 0.02 — 0.03 pc, 
A ~ 0.03 pc) objects that are located within cluster- forming cores in, for example , Ophi- 



uchus, Serpens, Perseus, and Orion (see discussion in lWard- Thompson et all 12007). Both 



the prestellar cores in e.g., Taurus, and the cluster-forming cores (which harbor multiple 
condensations) in the other clouds occupy only a very small fraction of the total volume 
and mass of their larger, more turbulent, molecular cloud complex fe.g. Jjohnstone et al 



20041 : iGoldsmith et all liooil ) . A compendium of current data from spectral line emis- 
sion, dust emission, and infrared absorption tends to show that cores exhibit central 

i ■ 1 1 I 1 I ; 1 I i 

density c oncentratio n (Ward- Thompson et al., 1994; Andre et al., 1996) subsonic i nward 
motions (jTafalla et all Il998t IWilliams et all Il999t \Lee et all l200lt ICaselli etail l2002h 
that often extend beyon d the nominal core boundary, near-un i form g as temperatures 
dBenson fc Mversl 19891). su bsonic internal turbulence (Myers, 1983 ; Fuller fe Mversl . 



1993c iGoodman et al.. 119981). and near-critical magnetic field strengths, when detected 



by the Zeeman effect ( Crutcher . 1999; Bo urke et al. . 2001 ) or inferred by the CF method 



(|Crutcher et all 120041) 



The measured subsonic infall motions constitute indirect evidence for a force that me- 
diates gravity. The mediative force may be due to the mag netic field , whos e strength in 
dense regions is very close to the critical value for collapse (jCrutcherl . Il999l ). If the frag- 
menting region has a subcritical mass-to-flux ratio, then the formation of cores will be 
regulated by ambipolar diffusion, i.e. neutral molecules diffusing past ions that are tied 
to magnetic fields, and will occur on a time scale that may significantly exceed the local 
dynamical time for typical ionization fractions. Conversely, if the dense gas is sup e rcriti- 
cal, then the fragmentation takes place on a dynamical time scale. iBasu fc Ciolekl (120041 
hereafter BC04) studied the nonlinear development of fragmentation instability in clouds 
that were initially either critical or decidedly supercritical. They found that supercritical 
fragmentation is characterized by somewhat supersonic motions on the core scale (~ 0.1 
pc) while critical fragmentation is characterized by subsonic inward motions on those 
scales. The nonlinear fragmentation of decidedly subcritical clouds is also presented in 
this paper as part of a comprehensive parameter study of the effects of mass-to-flux ra- 
tio, initial ionization fraction, and external pressure. We believe that the three broad 
categories of subcritical, transcritical, and supercritical fragmentation should all occur 
in the interstellar medium, and even within separate regions of a single molecular cloud 
complex. Our results can help to distinguish which mode of fragmentation is occurring 
in a given observed region. 

The linear theo ry of fragmentation o f a partially ionized, magnetic thin sheet has 
been presented bv lCiolek fc Basul (|2006l . hereafter CB06). An important quantity is the 
dimensionless critical mass-to-flux ratio, /i = 2nG 1 ^ 2 a/B, where a is the column density 
of a sheet and B is the strength of the magnetic field that is oriented perpendicular to 
the plane of the sheet. In the limit of flux- freezing, fragmentation can occur only if /x > 1, 
i.e. the sheet is supercritical. Conversely, gravitationally-driven fragmentation instability 
cannot occur at all if [i < 1, i.e. the sheet is subcritical. CB06 show that the inclusion of 
ambipolar drift in such a model means that fragmentation can occur for all mass-to-flux 
ratios, but on widely varying time scales and length scales. A very important result of 
CB06 is their Fig. 2 (see also Fig. 1 in this paper), which demonstrates that transcritical 
(/i w 1) fragmentation has a preferred scale that can be many times larger than Ax,m) the 
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wavelength of maximum growth rate in the thermal (i.e. nonmagnetic) limit. The latter 
is actually the preferred scale for both highly supercritical and highly subcritical clouds. 
This is because gravitational instability in subcritical clouds develops by ambipolar drift 
of neutrals past near-stationary magnetic field lines, and occ urs on the amb i polar diffusion 
time scale rather than the dynamical time scale (see also Langer , 19781 Mouschoviasl 



19781 ). The "resonant" transcritical modes with preferred scale A gjln 3> At.di occur due 



to a combination of ambipolar drift and field-line dragging; the latter leads to magnetic 
restoring forces that can stabilize the perturbations unless they are of the required large 
size. The time scale for ambipolar-diffusion mediated gravitational instability is also 
presented extensively in CB06. They showed that highly supercritical clouds undergo 
gravitational instability on the dynamical time i<j « Z /c s , where Z is the half-thickness 
(effectively the scale height) of the sheet and c s is the isothermal sound speed. On the 
other hand, highly subcritical clouds undergo instability on the quasistatic ambipolar 
diffusion time tad ~ 10 Z/c s for typical ionization fraction. Transcritical cl ouds undergo 
a hybrid instability on an intermediate time scale - see also IZweibell (|l998Jl for a similar 
result. Also, we note that the inclusion of no nlinear fluctuations can reduce the ambip olar 
diffusion time under certain circumstances ( Fatuzzo fc Adamsl . 2002 ; Zweibel . 20021 ); we 
do not model such effects in this paper. 

Our nonlinear simulations represent a significant extension of the parameter space of 
mo dels from that presented by B C04. It similarly extends the parameter space studied 
by llndebetouw fc Zweibel (2000), who presented nonaxisymmetric evolution of an in- 
finitesimally thin subcritical sheet, including the effects of magneti c tension but ignorin g 
magnetic pressure. Recent fully three-dimensional simulations by iKudoh et all (|2007h . 
including magnetic fields and ambipolar diffusion, have confirmed the basic results of 
BC04. However, that paper presents three representative models, and an extensive pa- 
rameter study remains computationally out of reach. Here, in addition to a parameter 
study, we carry out large numbers of simulations for each unique set of parameters. This 
is done in order to compile statistics on core spacings, mass distributions, and shapes. 
All simulations end at the time of the runaway central collapse of the first core, hence 
we are compiling information about the early phases of star formation in a molecular 
cloud. The subsequent history of the molecular cloud, after it has been stirred up by the 
initial star formation, remains to be determined. This paper represents the product of a 
total of over 700 separate simulations for models with 14 distinct sets of dimensionless 
parameters. 

An alternate approach to mode ling fragmentation is to input turbulence dire ctly into 
the dense thin-sheet model (e.g., Li fc Nakamural . l2004t iNakamura fc Lil . l2005h . rather 
than relegating its presence to an unmodeled envelope. This results in highly supersonic 
motions within the dense sheet itself, and more rapid formation of cores than in our 
models. We do not pursue that approach in the present study but leave it open to future 
investigation. 



1.2. Relation to Global Cloud Structure 



Our simulations represent an intermediate approach between attempting a global 
model of large-scale molecular cloud structure and of modeling the interior collapse of 
individual cores. Large-scale models need to account for the overall structure and self- 
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gravity of molecular clouds and cannot be realistically studied using a periodic-box model. 
Truly global three-dimensional models are computationally expensive and still rarely at- 
tempted. Smoothed particle hydrodynamics (SPH) techniques have been used su ccess- 
fully to model entire cluster-forming regions (jBate et al l l2003HBonnell et alllioolh . and 
can explain some important features of star formation like mass segregation, binary frac- 
tion, and the initial mass function. These models are not fully global in that they do not 
include the effect of the molecular cloud envelopes. They also do not include the effect 
of magnetic fields or feedback from outflows, so they canno t addr ess the observed low 



global star formation e fficiency (SFE) 1 — 5 % (jLada fc Ladal . I2003T ) . However, the latest 



version of such models (jPrice fc Batd . 120081 ) does include flux- frozen magnetic fields and 
yields somewhat lower SFE's than the non-magnetic models. Another three-dimensional 
approach which is fully global and based on finite-difference or SPH techniques is mod- 
eling the formation of a mol ecular cloud from large-scale supersonic gas flows and the 
early evolution of the c loud ( Vazquez-Semadeni et all 120061 2007 ; Heitsch et al. . 2008t 
Hennebelle et al. . l2008l ). These models yield an initially flattened cloud whose subsequent 



evolution is affected by several instabilities (thermal instability, thin-shell instability, and 
Kelvin- Hclmholtz instability). The overall structure evolves away from a sheet-like con- 
figuration, but individua l segments may be trea ted as such. Of the above mentioned 
studies, only the work of iHennebelle et af. ( 2008 ) includes the effect of magnetic fields, 
which may work to suppress some of the instabilities. It is not clear w hether this is an 
impor tant factor in that model. However, the smaller-scale simulation of iNakamura fc O 
(|2008l ) does show the formation of a sheet-like structure due to the presence of a dy- 
namically important magnetic field. Amongst our models, the cases with high bounding 
pressure may be most appropriate for the scenario of cloud formation due to colliding 
flows. Our low external pressure models are more appropriate to the alternate scenario 
where the flows that form a molecular cloud are driven by gravity and /or channeled by 
a large scale magnetic field, e . g. by the magneto- Jeans instability (|Kim et all 120021 ) or 
the Parker instability ( Parker . 19661 ). 

Despite the various theoretical emphases on flattened structure, we note that observa- 
tions of molecular clouds reveal complex morphologies with projected shapes that may 
not look like globally flattened structures. Molecular clouds have been described vari- 
ously in the lit erature as stratified ob jects supported by internally g enerated turbulence 
(McKee, 19991 ). or as fractal objects ( Elmegreen fc Falgaronel . Il996) in which the inter- 
nal pressure is not as relevant. Our local sheet model may be applicable to the dense 
subregions of the clouds (where stars actually form in weak or rich clusters) in either 
scenario . As an example from the first scenario describe d above, one-dimensional global 
models (|Kudoh fc Basil I2003L 120061 : iFolini et all l2004f ) of molecular clouds reveal that 
internally-driven turbulence yields large-amplitude motions in lower-density envelopes, 
while retaining transonic motions in embedded dense regions; the fragmentation of the 
latter may be described by our models. From an observational point of view, we may 
apply our models to dense star forming regio ns such as L1495 and HC1 2 in the Taurus 



molec ular cloud (see [Go ldsmith et al.l 



the L1688 cl uster-forming core in Ophi- 



uchus (iMotte et al.Lll998l ). or better yet to the Pipe Nebula (jMuench et alll2007t ). which 
represents an even earlier stage of evolution, with a large number of relatively quiescent 
prestellar cores that are found in a dense elongated region (the "stem" of the Pipe). This 
region may represent the best available laboratory for the study of the early stage of 
star formation, before feedback from star formation has significantly modified a cloud's 
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internal sructure and motions. There is also evidence that the stem of the Pipe N ebula 
is flattened along the direction of the mean magnetic field (jAlves fc Franco! . 120071) . This 
property is similar to the better established res ult for the elongated structures in Taurus 



(jGoodman et all ll99Ci iGoldsmith et all l2008h 



2. Physical Model 

We consider the evolution of weakly ionized, magnetic interstellar molecular clouds. 
The clouds are isothermal, having a temperature T. As presented in CB06, we model 
a cloud as a planar sheet or layer of infinite extent in the x- and y- directions of a 
Cartesian coordinate system (x, y, z). At each time t the sheet has a local vertical 
half-thickness Z(x,y,t). We take our model clouds to be thin: by this we mean that for 
any physical quantity f(x,y,z,t) the condition f/V p f 3> Z is always satisfied, where 
V p = xd/dx + yd/dy is the planar gradient operator. The magnetic field that threads a 
cloud has the form 

B z , cq (x, y, t)z for \z\ < Z(x, y, t), 

B(x,y,z,t) = <( B z (x,y,z,t)z (1) 
+B x (x, y, z, t)x + B y (x, y, z, t)y for \z\ > Z(x, y, i), 

where -B Z) eq is the vertical magnetic field strength in the equatorial plane. For \z\ — > oo, 
B — > B rc fZ, where B Te { is a uniform, constant reference magnetic field very far away 
from the sheet. The magnetic field components above the sheet can be determined from 
B z ,cq{x,y) at any time using the divergence-free nature of the magnetic field and the 
current-free approximation above the sheet (see CB06 for details). 

Some simplification is obtained by integrating the physical system of equations gov- 
erning the evolution (conservation of mass and momentum, Maxwell's equations, etc.) of 
a model cloud along the vertical axis from z = — Z(x, y) to z = +Z(x,y). In doing so, 
a "one-zone approximation" is used, in which the density and the x- and y- components 
of the neutral and ion velocities, as well as the x- and y- components of the gravita- 
tional field, are taken to be independent of height within the sheet. The volume density 
is calculated from the vertical pressure balance equation 

p Cs 2 = |G^ + P cxt + ^±^, (2) 

where P ex t is the external pressure on the sheet and B x and B y represent the values 
at the top surface of the sheet, z = +Z. This simplification is commonly referred to 
as the "thin-sheet approximation" ; the motivation for and the physical reasonability of 
it is discussed at length in Section 2 of CB06. It is the nonaxisymmetric extension of 
the axisymmetric thin-sheet models used to study ambipolar diff usion and gravitational 
collap se in magnetic interstellar c l ouds, as originally developed bv lCiolek fc MouschoviasI 
(|l993l ) and lBasu fc MouschoviasI (|l994T ). 
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2.1. Basic Equations 

We solve normalized versions of the magnetic thin-sheet equations as justified in CB06. 
The unit of velocity is taken to be c s , the column density unit is <7 n , , and the unit of accel- 
eration is 27rGer n .o, equal to the magnitude of vertical acceleration above the sheet. There- 
fore, the time unit is to = c s /2irGa n .o, and the length unit is Lq = cj:/2irGo- n fi. From this 
system we can also construct a unit of magnetic field strength, B — 2nG 1 ^ 2 c7 n fi- The 
unit of mass is Mq = Cg/(47r 2 G 2 <r ni o). Here, cr ni o is the uniform neutral column density 
of the background state, and G is the gravitational constant. With these normalizations, 
the equations used to determine the evolution of a model cloud are 



~~ oT" = "Vp ' (cr n V n ) , 

Ob 


(3) 


r\ 

-Q t {o-nVn) = - Vp • (a n v n v n ) + F T + F M + o- n g p , 


(4) 


— |p = ~ V p ' (Bz^qVi), 


(5) 


Ft = -CoffVpCTn, 


(6) 


F M = B z , eq (B p - Z WpB z , cq ) + 0(V p Z), 


(7) 


, T n iO ( Pn,o\ kl „ 
Vi=V n ^ '- — F M , 

O-n \ Pn J 


(8) 


r 2 2 (S^cxt + <JI) 

cff n (p^+«ir' 


(9) 




(10) 


z = ^-, 

2p n 


(11) 


9 P = -VpV, 


(12) 


V = ^" 1 [-^K)/fc z ], 


(13) 


Bp = -Vp*, 


(14) 


* = T- 1 [T(B z , cq -B Icl )/k z ] . 


(15) 



In the above equations, a n (x,y) = J_ z p n (x,y) dz is the column density of neutrals, 
B p (x,y) = B x (x 1 y)x + B y (x,y)y is the planar magnetic field at the top surface of 
the sheet, v n (x, y) — v x (x,y)x + v y (x,y)ij is the velocity of the neutrals in the plane, 
Vi(x, y) = Vi tX (x,y)x + Vi. y (x,y)y is the corresponding velocity of the ions, and the 
normalized initial mass density (in units of cr n ,o/^o) Pn.o = 4(1 + -Pcxt), where P oxt is 
defined below. The operator V p = x d/dx+y d/dy is the gradient in the planar directions 
within the sheet. The quantities ip(x,y) and &(x,y) are the scalar gravitational and 
magnetic potentials, respectively, also in the plane of the sheet. The vertical wavenumber 
k z = (fc^ + fc 2 ) 1 / 2 is a function of wavenumbers k x and k y in the plane of the sheet, and the 
operators T and T~ x represent the forward and inverse Fourier transforms, respectively, 
which we calculate numerically using an FFT technique. Terms of order 0(W P Z) in Fm, 
the magnetic force per unit area, are not written down for the sake of brevity, but are 
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included in the numerical code; their exact form is given in Sections 2.2 and 2.3 of CB06. 
All terms proportional to X? p Z are generally very small. 

We also note that the effect of nonzero V p Z and external pressure P ex t is accounted 
for in the vertically-integrated thermal pressure force per unit area, Ft, through the use 
of Cg ff . This can be seen by noting that 

+z +z 
F T = J V P P dz = V p J P dz- 2P cxt V p Z 
-z -z 

= 2V p (PZ-P cxt Z) , (16) 

where P is the pressure inside the sheet. In the above expression, we have used P = P cxt 
at the upper and lower surfaces of the sheet, and also that V p (+Z) = — V p (— Z). Using 
the ideal gas equation for an isothermal gas, P — /9 n c 2 , the expression for half-thickness 
(Eq. [11]), and the normalized equation for vertical hydrostatic equilibrium (Eq. [10] . 
where we ignore the relatively small term Bp for simplicity), it is straightforward to 
derive the normalized expression for Cg ff (Eq. [9]). 

The above equations contain the following dimensionless free parameters: P cxt = 
2-Pext I^Go\ o is the ratio of the external pressure acting on the sheet to the vertical 
self-gravitational stress of the reference state. The dimensionless neutral-ion collision 
time of the reference state, f n i.o = Tm,oAo; expresses the effect of ambipolar diffusion. In 
the limit f n i t o — > oo there is extremely poor neutral-ion collisional coupling, such that 
the ions and magnetic field have no effect on the neutrals. The opposite limit, f n i o = 0, 
corresponds to the neutrals being perfectly coupled to the ions due to frequent collisions, 
i.e. flux freezing. The neutral-ion collision time of the reference state is 

m\ + m„ 1 

r ni , = 1.4— ^ — — , (17) 

mi ni,o(aw)m 2 

where m\ is the ion mass, which we take to be 25 a.m.u., the mass of the typical 
atomic (Na + , Mg + ) and molecular (HCO + ) ion species in clouds, n^o is the ion num- 
ber density of the reference state, and (crw)iH, is the neutral-ion collisio n rate, equal to 
1.69 x 1(T 9 cm 3 s" 1 for H 2 -HCO + collisions (|McDaniel fc MasonL fl973h . The factor of 
1.4 in Eq. (jTTJ) accounts for the fact that the effect of helium is neglected in calculating 
the slowing-down time of the neutrals by collisions with ions. The parameter k\ is the 
exponent in the power-law expression that is used to calculate the ion density n- x as a 
function of neutral density n n , namely, 

ni - , (18) 

where we adopt fc, = 1/2 and IC (~ 10~ 5 cm ~ 3 / 2 ) for all models in this study (e.g., 



Elmegreen . 119791 : lUmebavashi fc Nakanol . [l980f ). but keep in mind that calculat ion of the 



ion chemistry network makes fc; a function of n n ( Ciolek fc Mouschoviasi 1998). Finally, 



£? rc f = B lc f/B = _B rc f /2nG 1 / 2 <T n fl is the dimensionless magnetic field strength of the 
reference state. For physical clarity, we use instead the dimensionless mass-to-flux ratio 
of the background reference state: 

Mo = 2^/^=4-1, (19) 



Br, 
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where (27rG 1 ^ 2 )^ 1 is the critical mass-to-flux ratio for gravitational collapse in our 
adopted thin-sheet geometry (CB06). Models with /j,q < 1 (B IC { > 1) are subcritical 
clouds, and those with fj,Q > 1 (-Brcf < 1) are supercritical. The initial mass-to- flux ratio 
is also related to the commonly-used plasma parameter 



ft= r£'°/« n + ft*)- (20) 

Typical values of our units are 



T 



1/2 



^ = ' 188 IitTkJ kms ' 1 ' (21) 

— *(l^f(^)- 

— -*(i^k)(^)- 

M o = 9 . 1 9x 1 O-3( Ij ^) 2 (l^)M , (24) 
B » = ««( I5 #i_) .G. (25) 

Here, we have used N nj o = cr n ,o/mn, where m n = 2.33 toh is the mean molecular mass of 
a neutral particle for an H2 gas with a 10% He abundance by number. Furthermore, we 
may calculate the number density of the background state as 

n n , = 2.31 x 10 5 (j^ps) 2 (l + PJ) cm-. (26) 

The dimensional background reference magnetic field strength for a given model is simply 
B Ic f = B /fi . Finally, the ionization fraction (= ri;/n n ) in the cloud may be expressed 
as 

= = 3.45 X 10- (£) (15!i=L!)" 2 ( 1 + P„)-" 2 . (27) 

2.2. Numerical Techniques, Boundary and Initial Conditions 

The system of Eqs. Q - (fT5|) are solved numerically in [x, y) coordinates using a 
multifluid non-ideal MHD code that was specifically developed for this purpose (BC04; 
CB06). Partial derivatives d/dx and d/dy are replaced with their finite-difference equiv- 
alents. Gradients are approximated using three-point central differences between mesh 
cells, while advect ion of mass and magnetic flux is prescribed by using the monotonic 



upwind scheme of Ivan Leerl (|1977h . Evolution of a model is carried out within a square 



computational domain of size L x L, spanning the region —L/2 < x < L/2 and —L/2 < 
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y < L/2. Typically, L is taken to be several times larger (up to a factor of 4) than the 
characteristic length scale of maximum gravitational instability A g . m (CB06; see, also, 
Section [3] below) . The computational domain is then divided into a set of iV 2 equally- 
sized mesh cells, each having an area L/N x L/N. Most of our simulations are run with 
L = 167r Lq and N = 128, and some have L = QAtt L$ and N = 512, so that the grid 
size Ax = Ay = 0.393 Lq in all cases. The mass resolution is then AM = 0.154 Mo, or 
1.42 x 1O~ 3 M using the standa rd values in Eq. l|2"i]). 

The numerical method of lines (|Schiesserl . ll99ll ) is employed, i.e. the first-order partial 
differential equations ((3j) - CSJ) are converted into a set of coupled ordinary differential 
equations (ODE's) in time, with one ODE for each physical variable at each cell. Hence, 
the system of ODE's has the form dy/dt = Q(y,t), where y and Q are both arrays 
of size VN 2 , V being the number of dependent variables. Time-integration of this sys- 
tem of OD E's is performed by using an Adams-Bashforth-Moulton predictor-corrector 
subroutine (|Shampind . [l99l . Numerical solution of Fourier transforms and inverse trans- 
forms, necessary to calculate the gravitational and magnetic potentials ip and \& at each 
time step (see Eq s. [13] and [15]), is done by using fast Fourier transform techniques 



(Press et al 



;e Kq s. 
19961) . 



Periodic conditions are applied to all physical variables at the boundary of the com- 
putational domain. The background reference state of a model cloud is characterized by 
a uniform column density <r n ,o and magnetic field B Zjeq fiZ = B rc f£. This means that the 
gravitational and magnetic forces are each identically zero in the uniform background 
state. The evolution of a model cloud is started at time t — by superposing a set of 
perturbations Sa n (x,y) that are random white noise with a root-mean-squared (rms) 
value that is 3% of cr n .o- To preserve the same local mass-to-flux ratio a n / 'B z ^ eq as in the 
uniform background state, initial magnetic field perturbations 6B zcq — (6cr n /<T n .o)B IC { 
are also introduced. 

Detailed tests of the accuracy of this MHD code were described in CB06. Full code 
runs were compared to exact linear solutions for the gravitationally unstable modes of 
thin-sheet magnetic clouds. The code was found to be in excellent agreement with these 
solutions. It correctly captured the temporal evolution of a model cloud in the linear 
regime of collapse, as exemplified by the growth time of the gravitational instability r g 
for a given fragmentation length scale A, for various values of the initial parameters /^o, 
f a ifi, and P e xt- In addition, we have run flux- freezing tests of subcritical models with 
random perturbations and verified that no spurious gravitational instability occurs in 
the absence of ambipolar diffusion. 

The principal motivation for our modeling clouds as thin sheets is that it significantly 
reduces the computational complexity of studying star formation, while still retaining 
many fundamental physical features necessary to understanding the dynamics of core 
formation and collapse within interstellar clouds. For instance, although model clouds 
are thin, they are not infinitesimally so, and we are able to incorporate both magnetic 
pressure and magnetic tension supporting forces (see Eq. [7]). Additionally, the vertical 
integration (along the direction of the z-axis) that is employed to derive the system 
of governing equations in the thin-sheet approximation (Eqs. [3] - [T5]) has the effect of 
turning the fully three-dimensional gravitational collapse problem into a computationally 
more tractable two-dimensional problem. As a result of these computational savings, our 
numerical code is able to run efficiently on a single workstation with minimal cpu times. 
Depending on the initial parameters, a full simulation can be completed in as little as an 
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hour or at most a single day. Hence, we are able to quickly generate a large number of 
models covering the entire physically relevant range of our free parameters, and produce 
a large quantity of models that can be used for statistical analysis. By contrast, three- 
dimensional MHD models require a dedicated workstation or a computer cluster, and 
their simulation completion times are orders of magnitude greater than that needed for 
our thin-sheet models. Our new code is written in the IDL programming language, which 
significantly speeds up the processes of code development, debugging, and visualization. 

We have also developed software to analyze the masses and shapes of cores arising from 
our simulations. In any snapshot of the evolution, we first isolate regions with column 
density (or mass-to-flux ratio in some cases) above some threshold value. Thresholds 
are usually chosen to be high enough that multiple peaks do not fall within a single 
contiguous region above the threshold. In cases where this happens, we have the option 
of manually isolating the cores. The mass of each core is found by adding the masses 
of each computational zone in the isolated region above the threshold. To determine 
the size and shape of a core we use the MPFITELLIPSE routine written in IDL by C. 
Markwardt, which returns the best-fit ellipse to the set of zones that constitute each 
core. The semimajor and semiminor axes of the best-fit ellipse, a and b, respectively, are 
obtained and used to determine the size s = yab and axis ratio b/a of each core. The 
vertical half-thicknesses Z of the zones within each core are averaged to find a mean 
half-thickness. The mean separation of fragments are found by locating, for each density 
peak associated with a core, the nearest peak of a neighboring core. These values are 
averaged over all cores in a simulation and over multiple model realizations. The periodic 
boundary conditions are also accounted for; we count any possible nearest neighbor that 
is just across the periodic boundary. 



3. Results 

3.1. Overview 

The efficiency of our two-dimensional code allows us to run a large number of simu- 
lations, with various combinations of the important parameters fj,o, T n ;.o, and P e xt- For 
each unique set of parameters we are also able to run a multitude of independent model 
realizations. Each realization is distinct in specific details, since the evolution is initiated 
by random (white noise) small-amplitude perturbations. However, the independent real- 
izations are statistically similar, and running a large number of models allows us to assess 
the level of randomness that contributes to distributions of various calculated quantities. 

Table 1 contains the parameters for each of ten models as well as the predicted mini- 
mum growth time T g , m , the associated wavelength of maximum instability A g , m , and the 
implied fragmentation mass M g . m = 7rcr nj oAg m /4, all obtained from the linear perturba- 
tion analysis of CB06. The first two quantities are in normalized form, but the masses 
have been converted to M & using the standard values in Eq. . Some insight into the 
numerical values of r g:m and A gim in Table 1 can be obtained by considering the growth 
rate and fragmentation scale of the fastest growing mode in the limit of no magnetic 
field, but finite external pressure. From the results of CB06, we find 
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r g (P cxt , B Iei = 0) = 2 ( 1 + 3 y 1/2 j£ = (i + ap^i/a ^ (28) 
(1 + P cxt ) c s c s 

A g (P cxt , B tei = 0) = ^ ^^ Lo = 2. (^-^ (29) 

(1 + Pextr V 1 + ^cxt / 



In the above equations, we have used the relation 
2L 

Zjn = = . 

(1 + Pext) 



(30) 



These results show that in the limit P cx t — + 0, the isothermal sheet has effective "Jeans 
length" Ax.m = — 2ttZq, and growth time Tr,m = 2Lo/c s = Zq/c s . The relation of 

these values to the sheet half-thickness and sound speed is intuitively understandable; the 
latter is essentially the dynamical time. The highly supercritical model 10, which also 
has low external pressure, does approach these limiting values of fragmentation scale 
and growth time. Equations ([28 ]) - ([29 ]) also bring out the interesting property that the 
fragmentation length and time scale of the isothermal sheet can be reduced significantly 
by increases in P ex t, but that the fragmentation scale co nverges to a fixed multiple of Zq 
(i.e. 6nZ ) as P ex t — » oo. This property was firs t noted by Elmegreen fc Elmegreen ( 19781 ) 



and studied further bv lLubow fc Pringld (|l993h . The significantly subcritical model 1 also 



has a fragmentation scale w 2ttZq, since instability occurs via neutral drift past near- 
stationary field lines; however the growth time r g , m w 10 Z /c s . The transcritical models 
3 and 4 have values of T Sjia more similar to the subcritical models than the supercritical 
ones, although their fragmentation scales A g , m are large. Models 9 and 10 have have 
Pcxt = 10, resulting in much smaller values of T gim and A gjm than corresponding models 
with P cxt =0.1, although the dynamically important magnetic field raises the values of 
both above the nonmagnetic limits. 

Our standard simulation box is four times larger in size than AT, m and more than twice 
Ag.m for most models. Exceptions to this are model 3 and model 7, and these simulations 
are carried out in larger simulation boxes of quadruple size (L — 6AttLq,N = 512) 
when collecting information on core properties. The total mass in the simulation box 
with L = 16ttL = 0.353 (iV n ,o/10 22 cm- 2 )- x (T/10 K) pc is M = 2.53 x 10 3 M = 
23.2 (iV ni o/10 22 cm" 2 )" 1 (r/10 K) 2 M Q . The expected fragmentation scales and masses 
for all models are significantly greater than the grid length resolution Ax and mass 
resolution AM quoted in Section 12.21 

Table 2 contains model parameters and key quantities at the end of each nonlinear 
model run. All runs end when cr n , max/c n .o = 10- This corresponds to a volume density 
enhancement p n /pn,o ~ 100 for models with P ext = 0.1 and p n /Pn0 ~ 10 for models 
with P cxt = 10. For each model, we list representative values of t rur G_I the physical time 
at which cr n , max = 10cr n ,0j l^nlmax, the maximum neutral speed in the simulated region 
at that time, and |i>i| max , the corresponding quantity for the ions. In general, models 
with relatively large values of fiQ, t^o, and P cxt tend to evolve faster than counterparts 



1 Experimentation with different rms amplitudes of the initial perturbation (so that Sa n /a n fi is in the 
range l%-6%), and variation of the power spectrum away from white noise but with the fixed standard 
rms value, reveal that the values of t run can vary in the range 10%-20% from the values quoted in Tabic 
2. 



12 



Table 1 

Summary of Parameters and Results of Linear Theory 



Model 


MO 




■Pcxt 


T g,m 


<^g,m 


M g,m 


1 





.5 


0.2 


0.1 


20.2 


14.3 


1.48 


2 





.8 


0.2 


0.1 


17.6 


16.5 


1.96 


3 


1 


.0 


0.2 


0.1 


14.3 


24.3 


4.26 


4 


1 


.1 


0.2 


0.1 


10.8 


54.9 


21.8 


5 


2. 


.0 


0.2 


0.1 


3.21 


23.9 


4.12 


6 


10. 


.0 


0.2 


0.1 


2.11 


13.9 


1.40 


7 


1. 


.0 


0.1 


0.1 


27.5 


28.5 


5.86 


8 


1. 


.0 


0.1 


0.1 


7.93 


20.0 


2.89 


9 





.5 


0.2 


10.0 


4.87 


3.43 


0.09 


10 


1 


.0 


0.2 


10.0 


3.33 


4.99 


0.18 



Times and lengths are normalized to to and Lo, respectively. Masses are converted to Mq using Eq. 

with smaller values of these parameters. For the gravity-dominated models (P C xt = 0.1), 
increasing values of fio and/or f n ;.o result in greater values of |f n |max- The values of 
|^i|max illustrate that the systematic motions within the nonlinearly developed cores are 
gravitationally driven, so that ions lag behind neutrals somewhat. However, the difference 
between the speeds of the two species remains less than 0.1 c s . We perform an analysis 
of core properties as described in Section 12.21 after defining a core as an enclosed region 
with a n /o~ n ,o > 2 at the end of the simulation. Since each simulation typically yields only 
a handful of cores, the models are run a large number of times to generate significant core 
statistics. Models 1, 2, 3, and 5 were run 100 times each, while models 4 and 7, which 
need to be run on an expanded grid due to very large fragmentation scales, were run 22 
and 6 times, respectively. Model 8 ran 50 times and models 9 and 10 were run 25 times 
each. Several averaged properties of the resulting cores are presented in columns 8-12 of 
Table 2. These quantities are (A), the average distance between cores, (M), the average 
mass within a core (converted to M using Eq. [24]), (s), the average size of a core, (b/a), 
the average axis ratio of a core, and (Z), the average value of the half-thickness of a core. 

Fig. [U shows the preferred fragmentation scales A gjm versus fio from the linear analysis 
of CB06 for two separate values of the dimensionless external pressure (P cx t = 0.1 on top 
and P cxt = 10 below). Both lines represent models with f n ;.o = 0.2, which is our adopted 
standard value based on typical observationally-inferred ionization levels (see Eq. [27] 
above and Eq. [29] of CB06). Overlaid on each solid line are average core spacings (A) 
(see Table 2) calculated from the nonlinear endpoint of our simulations. Each data point 
represents an average of core spacings for a large number (20 to 100) of simulations. 
The lower line contains (A) data for up to 25 runs each of an additional four models 
with P cxt = 10 that do not, for the sake of brevity, have their data compiled in Tables 1 
and 2. These results show that the linear theory can be used with confidence to predict the 
average fragmentation properties of clouds even in a fully nonlinear stage of development. 
The tabulated fragmentation scales are slightly below the predictions of linear theory in 
the range fio ~ 1 — 2, for P ext = 0.1. This is due to occasional subfragmentation of the 
initially large (irregularly shaped) fragments as they become decidedly supercritical. We 
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Fig. 1. Preferred fragmentation scales \ Sl m from linear stability analysis llCio lek & Basu, 2006) compared 
with the average spacing of density peaks in the nonlinear simulations. The upper solid line is the 
calculated dependence of the wavelength with maximum growth rate A gjm versus no for fixed parameters 
f n i = 0.2 and P ex t = 0.1. The lower solid line is the same but for P cx t = 10. The triangles represent the 
average spacing of density peaks (with peak cr n > 2cr n! o) tabulated from a large number of simulations 
at each of fio = 0.5, 0.8, 1.0, 1.1, 2.0, and 10.0 with f n i o = 0.2 and P cx t = 0.1. The squares represent the 
same but with P cx t = 10. 
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Table 2 

Main Results for Model Clouds 



Model 


Mo 


7ni,0 


-Pcxt 


trun 


| | max 


|f i |max 


w 


(M) 


<s) 


(b/a) 


(Z) 


1 


0.5 


0.2 


0.1 


204 


0.39 


0.35 


15.0 


1.26 


1.25 


0.74 


0.71 


2 


0.8 


0.2 


0.1 


167 


0.62 


0.55 


17.7 


1.56 


1.49 


0.81 


0.72 


3 


1.0 


0.2 


0.1 


121 


0.70 


0.64 


20.1 


3.26 


2.09 


0.69 


0.65 


4 


1.1 


0.2 


0.1 


88 


0.95 


0.90 


47.1 


6.62 


3.23 


0.66 


0.75 


5 


2.0 


0.2 


0.1 


23 


1.1 


1.0 


19.1 


3.41 


2.08 


0.57 


0.67 


6 


10.0 


0.2 


0.1 


12 


1.2 


1.2 


12.8 


0.99 


1.04 


0.53 


0.70 


7 


1.0 


0.1 


0.1 


261 


0.66 


0.62 


31.4 


3.29 


2.13 


0.77 


0.78 


8 


1.0 


0.1 


0.1 


61 


0.73 


0.63 


18.5 


2.02 


1.69 


0.67 


0.71 


9 


0.5 


0.2 


10.0 


44 


0.70 


0.63 


4.0 


0.34 


0.60 


0.62 


0.072 


10 


1.0 


0.2 


10.0 


22 


0.50 


0.38 


5.3 


0.46 


0.74 


0.60 


0.074 



Times and lengths are normalized to to and Lo, respectively. Speeds are normalized to c s . Masses 
are converted to Mq using Eq. 1241 ). Core data for models 4 and 7 are compiled from runs with 
N = 512, L = 64ttL . 

return to this issue when discussing the models with fiQ = 1.1. 



3.2. The Effect of Varying /iq 



3.2.1. Time Evolution 

Fig. [2] shows the time evolution of the maximum neutral column density o~ n max and 
maximum mass-to-flux ratio /^ ma x, both normalized to their initial values, for models 1, 
3, and 5, which have [io = 0.5, 1.0, and 2.0, respectively, and fixed parameters f n i,o = 0.2 
and P cxt — 0.1. These three models have r gim /io = 20.2, 14.3, and 3.2, respectively. The 
actual time t run to runaway collapse of a core, starting from small-amplitude white-noise 
perturbations, is about 7—10 times T g m for these models. Review of Table 1 and Table 
2 shows that i r un/ T g,m ~ 7 — 10 is a generic feature of all our models. The clouds with 
initial critical and subcritical mass-to-flux ratio have a prolonged period of dormancy, 
compared to the supercritical models. This is due to the need for ambipolar diffusion 
to operate before collapse sets in for both cases. The dashed lines in Fig. [2] show how 
much the mass-to-flux ratio changes during the evolution. The initially subcritical cloud 
requires a significant increase of /i m ax before collapse begins. Although /i max appears to 
be diverging at the end of the simulations, it is in fact increasing much more slowly than 
On. max, and will not asymptot ically diverge. This can be se en in previously published 
work i.e. in Fig. 2 of each of Ciolek fc Mouschoviasl ( 1994 ) and Basu fc Mouschoviasl 
(|l994h . 



3.2.2. Column Density and Velocity Structure 

Fig. E] shows the column density map and velocity vectors of neutrals at the end of 
the simulations for parameters /io = 0.5,0.8, 1.0, 1.1,2.0, and 10.0, with f n io = 0.2 and 
fcxt = 0.1. All simulations end when crn.max/cn.o = 10, but the time at which this is 
reached is different in each model. These times for the various models (in order of increas- 
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Fig. 2. Time evolution of maximum values of surface density and mass-to-flux ratio in three simulations. 
The solid lines show the evolution of the maximum value of surface density in the simulation, o" n ,max/o"n,o, 
versus time t/to- This is shown for models 1, 3, and 5, which have initial mass-to-flux ratio values /to = 
0.5, 1, and 2, respectively. For each model, a dashed line shows the evolution of the maximum mass-to-flux 
ratio in the simulation, /ti ma x, normalized to /no. 

ing fi ) are t/t = 203.7, 166.5, 120.9, 88.1, 22.8, and 12.4. There is a striking variation in 
the spacings of cores as [i changes. The highly subcritical case fi — 0-5 fragments on 
essentially the nonmagnetic preferred scale A^m = 47rL , since the evolution is charac- 
terized by diffusion of neutrals past near-stationary magnetic field lines, i.e. a Jeans-like 
instability but on a diffusive time scale. As predicted by the linear theory (CB06), there 
is a peak in the fragmentation spacing near = 1- The peak occurs at /i = 1.1 when 
Tni.o — 0-2 and P oxt = 0.1. The predicted fragmentation scale is A gim = 4.2 At, m , and 
indeed our simulation with box width 4 Ax, m yields only one core, although it seems to 
be undergoing a secondary fragmentation into two pieces. 

The velocity vectors of the neutral flow arc normalized to the same scale in each 
frame, and the horizontal or vertical spacing of the footpoints is equal to 0.5 c s . There is 
a monotonic increase of the typical neutral speeds as fi increases (see values of |w n |max 
in Table 2). The supercritical models, unlike their subcritical and critical counterparts, 
have large-scale flow patterns with velocities in the approximate range (0.5 — 1.0) Cg, 
and maximum speeds associated with the most fully developed cores that are mildly 
supersonic at distances ~ 0.1 pc from the core centers. Interestingly, the essentially 
hydrodynamic model with fi — 10 has only somewhat greater systematic speeds than 
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the model with /irj = 2. This is because the fragmentation scale is smaller, so that each 
core has a weaker gravitational influence on its surroundings at this stage of development. 
However, recall that the frames are at different physical times since the models with 
smaller values of reach cr Uimax /cr nt o = 10 at progressively later times. Spatial profiles 
of v n in the vicinity of cores are presented in BC04, for some models, and we do not 
present them again in this paper. The trend of values of maximum neutral speed |v n | max 
and maximum ion speed |i>i| ma x (Table 2) reveal the predictions of our models for the 
observable motions on the core scale. The infall motions are gravitationally driven, so 
that the ion speed lags the neutral speed in all cases. The overall rms speed in the 
entire simulation region is always quite small for both neutrals and ions. The rms speed 
of neutrals, w n ,nns> falls in the range (0.05 — 0.21) c s for the various models, and the 
corresponding quantity v^rms in the range (0.04 — 0.20) c s . 

Since the standard box size L = IQttLq allows only one fragment to form initially in 
the model with fia = 1.1 (as seen in Fig. [3]), we ran another model with four times larger 
box size but the same resolution, so that L = 6AttLq and N = 512. This allows the 
formation of multiple fragments, since the preferred fragmentation scale from the linear 
theory is A gim = 54.9Lo- Fig.rj]shows the column density and velocity vectors at the end 
of one such simulation. Velocity vectors are again normalized such that the horizontal or 
vertical spacing between footpoints equals 0.5 c s . This larger simulation does show that 
multiple fragments form with spacings approximately as predicted by the linear theory, 
but that there is also a tendency for cores to subfragment into two density peaks. An 
analysis of the result of 25 separate simulations with different random realizations of the 
initial state reveals that the average distance between density peaks is 47.1Lo, which is 
somewhat smaller than A g!m = 54.9Lrj. This is explained by the occasional presence of 
secondary density peaks within the initially formed fragments. While most model clouds 
undergo single-stage fragmentation into essentially thermal critical (Jeans-like) fragments 
of size ~ Ax,m, the transcritical clouds form first-stage fragments many times larger than 
AT,m) followed by a possible second-stage fragmentation when the mass-to-flux ratio of 
the fragment becomes decidedly supercritical due to ambipolar diffusion. This second- 
stage fragmentation may be favored because both the preferred fragmentation scale and 
growth time of gravitational instability drop precipitously as a cloud makes the transition 
from transcritical to supercritical (see Figs. Id and 2 of CB06). The initially rather large 
fragment may itself be prone to fragmentation because of its irregular shape. 

Fig. [U shows an alternate view of the column density, using surface plots at the end 
of simulation runs, with parameters corresponding to those of models 1, 3 and 5. These 
models start from different realizations of the initial state than those of the models 
presented in Fig.[3l Density peaks occur in different locations but represent an equivalent 
outcome statistically. Animations of the time evolution of the surface plots are available 
onlinj^l 



2 The animations of the models shown in Figs. [5] and [6] reveal that they reach the final state with 
"n/fn.o = 10 in a time 15-20% less than that quoted in Table 2. This is because the initial perturbation is 
white noise but with the smallest wavelengths (A < 4 grid cells) damped out. This results in slightly more 
power in the longer wavelengths (for a fixed rms perturbation level) , including the preferred fragmentation 
scale A gl m, and a consequent quicker development of the favored mode. 
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Fig. 3. Image and contours of column density a n (x, j/)/o" nj o, and velocity vectors of neutrals, for six 
different models at the time that <T njmax /(j n> o = 10. All models have f n ; q = 0.2 and P e xt = 0.1. 
Top left: fi = 0.5. Top right: fj, = 0.8. Middle left: fi = 1.0. Middle right: fi = 1.1. Bottom left: 
fj-0 = 2.0. Bottom right: [iq = 10.0. The color table is applied to the logarithm of column density and 
the contour lines represent values of a n /cr n fi spaced in multiplicative increments of 2 1 / 2 , having the 
values [0.7, 1.0, 1.4,2, 2.8,4.0,...]. The horizontal or vertical distance between footpoints of velocity vectors 
corresponds to a speed 0.5c s . We use the normalized spatial coordinates x' = x/Ax, m and y' = i//\t mi 
where At m is the wavelength of maximum growth rate in the nonmagnetic limit with P cx t = 0. 
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Fig. 4. Column density and velocity vectors as in Fig. [3] but with model 4 parameters (fio = 1.1) that 
is run with four times larger computational region in each direction. The simulation region consists of 
512 X 512 zones. The horizontal or vertical distance between footpoints of velocity vectors corresponds 
to a speed 0.5 c s . 




Fig. 5. Surface plot of column density a n (x,y)/a n: o at the end of the simulation for models 1, 3, and 
5 with (left to right) fio = 0.5,1.0,2.0. Animations of the time evolution for each model are available 
online. 
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3.2.3. Magnetic Field Lines 

For the force- free and current-free region above our model sheet, the three-dimensional 
structure of the magnetic field is obtained by solving Laplace's equation for the scalar 
magnetic potential ^m(x, y, z). The two-dimensional Fourier Transform of the magnetic 
potential at any height z above the sheet is related to its known planar value ^(x,y) = 
*m(x, y,z = 0) via 

* M (x, y, z) = T- 1 {F[*{x, y)} exp(-fc z \z\)} . (31) 

In the above expression, T {J-^ 1 ) represents the forward (inverse) Fourier Transform 
in a two dimensional (x, y) plane for a fixed z, ^ is obtained from Eq. (I15p . and k z = 
{k 2 + k 2 ) 1 / 2 . We use the FFT technique to efficiently calculate at any height z, 
which then leads to the various components of the magnetic field above the sheet via 
B — B le {Z = — V^m (see CB06 for a justification of this expression). 

Fig. [6] shows images of the final state column density for two models with /xq = 1.0 
and 2.0, respectively. These are also independent realizations and differ in specific details 
from models presented earlier. Overlaid on the column density images are the magnetic 
field lines extending above the sheet. These are three-dimensional images of a sheet plus 
field lines above, viewed from an angle of approximately 10° from the direction of the 
background magnetic field. Clearly, the supercritical model has more curvature in the 
field lines, as the contraction proceeds primarily with field-line dragging. The critical 
model produces its cores via a hybrid mode including both neutral-ion slip (ambipolar 
diffusion) and field-line dragging, hence the the lesser amount of field line curvature. The 
relative amounts of field line curvature in the cloud and within dense cores are quantified 
by calculating the quantity 9 = tnn -1 (\B p \/ B z<eq ), where \B p \ = (B 2 + B 2 ) 1 / 2 is the 
magnitude of the planar magnetic field at any location on the sheet-like cloud. Hence, 
9 is the angle that a field line makes with the vertical direction at any location at the 
top or bottom surface of the sheet. To quantify the differences in field line bending from 
subcritical to transcritical to supercritical clouds, we note that models 1, 3, and 5, with 
/io = (0.5,1.0,2.0), have average values # av = (1.7°, 8.3°, 18°), and maximum values 
(probing the most evolved core in each simulation) # max = (20°, 30°, 46°). For the model 
with [io = 0.5, 9 max is comparable to that presented in Fig. 2 of lCiolek ( 19961 ). 



Our ability to model fragmentation with varying levels of magnetic support and neutral- 
ion slip opens up the possibility of making a detailed com parison of observ ed hourglass 
morphologies of magnetic field lines, where measured (e.g. Schleuninel . 19981 ). with theo- 



retical models so that the curvature of field lines may be used as a proxy to measure the 
ambient magnetic field strength. 



3.2.4. Shapes 

The core identification techniques described in Section 12.21 yield core sizes, shapes, 
and masses, whose average values are tabulated in Table 2. Here in Fig. [7| we present 
a more detailed histogram of core shape distributions for each of models 1 through 6. 
This isolates the effect of mass-to-flux ratio on core shapes. The statistics are generated 
by running each of model 1-3 and 5-6 for 100 distinct realizations with the usual box 
size (L = 16ttLq,N = 128). Model 4 is run 22 times due to the larger box size [L = 
64wLq,N = 512) necessitated by the large fragmentation scale. The number of cores 
generated are 547, 367, 187, 126, 272, and 659, for models 1 through 6, respectively. The 



20 



Fig. 6. Image of gas column density cr u (x, j/)/<r n ,o and superposed magnetic field lines for models 3 and 
5, with fio = 1.0 (left) and fio = 2.0 (right). The magnetic field lines extend above the sheet, and the 
image and field lines are seen from a viewing angle of about 10°. Animations of the evolution of the 
column density are available online. The field lines appear in the last frame of the animation. 



different numbers reflect the varying numbers of cores arising for each set of parameters 
(see Fig. [3]) as well as the smaller number of runs for model 4. However, in each case we 
have sufficient numbers to make inferences about the differences between models. 

Fig. [7] reveals in detail what is apparent by a visual inspection of Fig. [3] The core shape 
distribution contains many circular objects (axis ratio b/a ~ 1) for subcritical clouds, 
but contains progressively more elongated cores as /io increases. All initially supercritical 
models have a peak axis ratio that is distinctly non-circular. All objects are also flattened 
in the vertical direction and usually have the shortest dimension along the magnetic field 
(see values of (Z) in Table 2). The underlying physical explanation is that quasistatic 
formation of cores (for /io < 1) allows for growth in all directions more equally, whereas 
dynamical gravity-dominated formati on will strongly accent uate the anisotropics (k x ^ 



by) that are present at the outset (see lMivama et al 



The supercritical models have a mean axis ratio in the sheet plane ~ 0. 5, which is in 



rough agreement with the observed mean projected axis ratio of dense cores (jMvers et al 



1991). However, a deprojection of t he observed axis ratios yields intrinsic three-dimensiona l 



shapes that are inh erently triaxial (jJones et al.l . l200lUJones fc Basull2002tlGoodwin et al 



2002 ; Tassisl 2007 ). with mean axis ratios b/a~ 0.9 and c/a rj 0.4 — 0.5. Since the direc- 
tion of the smallest axis (c) corresponds to our preferred direction of flattening (z), the 
deprojected b/a values can be compared directly with our models. We find reasonable 
agreement for the subcritical models 1 and 2. The deprojected values of c/a can also be 
compared with our thin-sheet models, in which effectively c/a — (Z) /a — (Z) y {b/a) /(s) . 
There is reasonable agreement here again for the subcritical models 1 and 2, as well as 
for the highly supercritical model 6, which have, respectively, c/a — 0.49,0.43, and 0.49. 



3.2.5. Core Mass Distributions 

In Fig.[8]we present histograms of core mass distributions generated from the multiple 
runs of models 1 through 6, as described in Section 13.2.41 Each core is defined as an 
enclosed region with o~ n /o~ n ,o > 2 that is present at the end of the simulation, when 



21 




0.15 0.55 0.95 0.15 0.55 0.95 0.15 0.55 0.95 

Axis Ratio Axis Ratio Axis Ratio 

Fig. 7. Histograms of axis ratios b/a of best fit ellipses to dense regions with ir n /Vn,0 > 2, measured at 
the end of simulations with fio = 0.5, 0.8, 1.0.1.1, 2.0, 10.0 as labeled, corresponding to models 1 through 
6 in Table 2. Each figure is the result of a compilation of results of a large number of simulations. The 
bin width is 0.1. 

<T n,max/ <T n,o = 10- For comparison with observations, we have converted our calculated 
masses to Mq using an assumed background number column density N n ,o = 10 22 cm -2 
and temperature T = 10 K (see Eq. [24]). 

We note that the variation of the peak masses (and the average masses tabulated in 
Table 2) from one model to another are in qualitative agreement with the predictions of 
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linear theory (Table 1). Furthermore, the /io = 0.5 and fig = 10 models generate very 
similar core mass distributions that are difficult to distinguish. This is not surprising 
since they have such similar preferred fragmentation scales. 

The striking feature of each of the histograms is the very sharp descent at masses 
greater than the peak of the distribution. Gravitational fragmentation yields a very strong 
preferred mass scale. The peak value itself is more ambiguous and can vary according to 
the magnetic field strength, the background column density, the cloud temperature, and 
the contour level we use to define the core. In contrast, the slope on the low-mass side is 
much shallower. This is due to the capture of emerging cores at the end of any simulation. 
Many of those cores are expected to grow in time and move over to the right by the time 
their peaks undergo runaway collapse and form a star. The distribution may be described 
as relatively narrow and lognormal-like, but with a broader tail at the low-mass side due 
to the temporal spread of core ages. 

The steep decline of the of the mass distributions beyond the peak (d log N/ d log M w 
—5 is typical) in our study is in c ontra st to that observed for condensations in cluster 
forming regions (e.g. Motte et al. . 1998f ). where d log N/d log M ~ —1.5 at high masses. 
Gravitational fragmentation under the conditions studied in this paper and at the time 
snapshot chosen here yields a very strong preference for a characteristic mass. We discuss 
possible mechanisms of broadening the mass distribution in Section [4j 



3.2.6. Supercritical Cores 

For clouds that start with subcritical or transcritical initial conditions, there is avail- 
able a more physical definition of a "core" , i.e. a region that is significantly supercritical 
and enclosed within a subcritical common cloud envelope. Axisymmetric simulations 
of cores that evolve initially by ambipolar drift have shown that the contraction be- 
comes very rapid by the time that fj, sa 2 in the central region , leaving behind a more 
slowly evolving and essentially subcritical envelope (see, e.g . Fiedler k. Mouschoviasl 
Il993t ICiolek fc Mouschoviaslll994 iBasu fc Mouschoviasl[l99l . 

Fig.[9]shows images and contour maps of fi(x, y) at the end of the simulation for models 
1 and 3, respectively. The initially subcritical (/io — 0.5) model 1 has peaks in n(x,y) 
coinciding with the major peaks in a(x,y) (see Fig. [3] upper left). However, note that 
the density condensations are either largely or even entirely subcritical (/i < 1) at this 
stage. This image shows that subcritical clouds can have observable density enhancements 
which may still be partially or entirely subcritical, because they are still in the process of 
ambipolar-drift-driven gravitational instability. The image and contours for the initially 
critical (/io = 1.0) cloud shows that the cloud naturally separates into supercritical and 
subcritical regions, due to ambipolar diffusion and a fixed total magnetic flux threading 
the cloud. The newly created supercritical regions are extended and typically contain 
more than one density and mass-to-flux ratio peak within them. In this case, all density 
peaks are associated with gas that has fi > 1. 

The presence of supercritical regions embedded within a common subcritical envelope 
allows us to define cores in a more physical way than the previous definition as regions 
with o~ n /<J n fl > 2. The latter is a somewhat arbitrary designation, as indeed are all 
observational definitions of cores. However, the definition of supercritical cores has its 
own ambiguities, as a simple definition of regions with /i > 1 yields extended regions in 
model 3 with multiple density peaks. We find that a viable working definition is that a 
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Fig. 8. Histograms of masses contained within regions with cr n /cr ni o > 2, measured at the end of simula- 
tions with fio = 0.5,0.8, 1.0.1.1,2.0, 10.0 as labeled. Each figure is the result of a compilation of results 
of a large number of simulations. The bin width is 0.1. 



core is a region with /i > 1.3. This isolates individual density peaks in both models, and 
is consistent with the earlier axisymmetric findings that a mass-to-flux ratio somewhat 
above the critical value is necessa ry before rapid collapse and sep aration from the envelope 



above tne critical value is necessa ry betore rapid collapse and s eparation from trie envelope 
becomes apparent. For example, Ciolek fc Mouschoviasl ( 1993f ) found that /i > 1.23 was 
required for the absence of any available axisymmetric equilibrium state, using a similar 
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Fig. 9. Image and contours of /i(x,y), the mass-to-flux ratio in units of the critical value for collapse. 
Regions with /i > 1 are displayed with a color table, while regions with /i < 1 are black. The contour 
lines are spaced in additive increments of 0.1. Left: Final snapshot of simulation with fig = 0.5. Right: 
Final snapshot of simulation with no = 1.0. 

value of P cx t as we do. We compile data from 100 runs of each model with distinct random 
realizations of the initial states and present the core mass distribution for each of models 
1 and 3 in Fig. 1101 The conversion to dimensional masses is done in the same manner as 
for Fig. [5] The resulting distributions have a peak mass that is somewhat smaller than 
found using the different core definition used for Fig.[5J However, these distributions also 
have a very sharp decline at higher masses. 

3.3. The Effect of Varying f n i t o 

Eq. (|27|) shows that the ionization fraction x\ at a given neutral density n n increases 
(decreases) linearly as f^Q decreases (increases). We investigate the effect of decreasing 
and increasing f^o by a factor of two from its standard value (which can be accomplished 
by changing the factor JC in Eq. [TB]) in models 7 and 8, respectively. The other two 
parameters are kept fixed at fio — 1.0 and P cx t =0.1. The characteristic growth times 
of instability T giin scale approximately cx t^q oc x\$ (see also Table 1), where x- x $ is the 
ionization fraction at the background number density n n ,o- Models 7, 3, and 8 have f n i,o = 
(0.1, 0.2, 0.4), r g:in = (27.5, 14.3, 7.9) x t , and A g , m = (28.5, 24.3, 20.0) x L , respectively. 
Furthermore, Table 2 reveals that the time to runaway collapse (cr njmax /(7 ni o > 10) is 
~ 10T g m when starting from small-amplitude white noise perturbations, as generally 
found in our parameter study. This means that our high ionization-fraction model 7 has 
the largest value of t Iun (= 261 to) m our parameter study. This also leads to the largest 
age spread of cores in any of our simulations. This is measured by the fact that a typical 
simulation, when run (for compiling statistics) in a large box with L = 6Att,N = 512, 
has many cores that are just beginning to emerge when the first core goes into a runaway 
collapse. Hence, our value for (A) is calculated with a lower core threshold cr n /<7 n ,o > V2 
for this one model. Our analysis reveals that the fragmentation scales in the nonlinear 
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Fig. 10. Histogram of masses contained within regions that are significantly supercritical, specifically 
fi > 1.3, measured at the end of simulations with fj,o = 0.5 and fio = 1.0. Each figure is the result of a 
compilation of results of a large number of simulations. The bin width is 0.1. 

phase are indeed comparable amongst models 3, 7, and 8, and in good agreement with 
the linear theory prediction. We conclude that the effect of varying ionization (for a fixed 
/^o) within the range studied is primarily in the rate of evolution. 

Fig.HUshows images and contours of the density, as well as velocity vectors, for models 
7 and 8. The time to reach runaway collapse is about four times longer for model 7 than 
for model 8, consistent with its value of f n i o being four times smaller. Model 8 has 
Tni.o = 0.4, hence poorer neutral-ion coupling and therefore reduced magnetic support. 
This results in slightly greater infall speeds, slightly larger number of fragments, and 
cores which are slightly more elongated. These are all consistent with its more dynamical 
evolution. 

3.4. The Effect of Varying P ext 

Our models with P oxt = 10 may represent the effect of pressured environments such as 
sheets brought together by the presence of shocked gas (e.g. stellar winds or supernovae) 
and being embedded in or adjoining an H II region. These models could represent an 
example of "induced" star formation in a manner related but not equivalent to that of 
an initial turbulent flow with high ram pressure. 
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Fig. 11. Column density and velocity vectors as in Fig. [3] but for models with T n i o = 0.1 (left) and 
f n ; o = 0.4 (right). Both models have fio = 1.0 and P cx t = 0.1. 



Models 9 and 10 both develop extreme clustering in comparison to the other models. 
Table 2 shows that the fragmentation scales are about 1/3 to 1/4 of that for the corre- 
sponding models with the same mass-to- flux ratio but small P xt- The fragments grow 
initially through a pressure-driven mode and the spacing is in excellent agreement with 
the predictions of linear theory (Fig. [T] and Table 1). However, our results show that the 
nonlinear instability does develop into a gravitationally-driven runaway collapse. The 
maximum speeds are still subsonic at the end of our simulation, due to the relatively 
weak gravitational influence of each compact core. As well as having the smallest frag- 
mentation scales, these models also have the shortest time scales to runaway collapse. 
For the fiducial N n = 10 22 cm' 2 , models 9 (/Uq = 0.5) and 10 (fio = 1.0) have values 
of < mn = (1.6 Myr,0.81 Myr) and (A) = (5800 AU, 7700 AU), respectively. Both sets of 
numbers are considerably smaller than for the models 1 and 3, which have corresponding 
values of hq but P e xt =0.1. Fig. [T2l shows the clustering properties of model 10 at the 
end of the simulation, which is very similar to the corresponding image for model 9 (not 
shown). This fragmentation model clearly produces a much richer cluster than in the 
relatively unpressured environments presented earlier. Velocity vectors are not shown in 
this image due to confusion arising from infall onto so many peaks. A careful inspec- 
tion of the image reveals a variety of core spacings and sizes at this stage of evolution. 
Interestingly, the average core spacing (A) = 5.3L is in very good agreement with the 
preferred wavelength in linear theory, A gjm = 5.0-Lo- These length scales are well resolved 
in our simulations. However, the average core mass (M) significantly exceeds the linear 
theory value Af g m . Only models 9 and 10 show such a large discrepancy between these 
values. We attribute it to the very small sizes of the cores (see (s) values in Table 2) 
in these simulations. This means that the cores themselves are barely resolved and the 
mass estimates should be taken as approximate values that likely represent upper limits. 
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Fig. 12. Column density as in Fig. [3] but for a model with Poxt = 10. Other parameters of this model 
are fio = 1.0 and T n i o = 0.2. 

4. Discussion 



Our simulations of the nonlinear development of gravitational instability under the 
influence of magnetic fields and ambipolar diffusion start from a background state of 
uniform column density and magnetic field strength. Small-amplitude white-noise per- 
turbations initiate the evolution and eventually lead to the nonlinear growth of fragments. 
Averaging over a large number of simulations reveals that the average spacing of non- 
linearly developed cores is essentially that predicted from the preferred fragmentation 
scales in linear perturbation theory (CB06). However, the time to reach fully developed 
runaway collapse is up to ten times longer than that of the eigenmode with minimum 
growth time T gjm . The quantity r gjm itself varies from sa Zq/c s (essentially the free-fall 
time « 1/ y/Gpnfl for unpressured sheets) for highly supercritical models to w 1QZq/c s 
for highly subcritical models (for a typical neutral- ion coupling level). The times to reach 
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runaway collapse vary widely amongst models with different mass-to-flux ratios, ioniza- 
tion fractions, and external pressures. For a cloud with N D £ = 10 22 cm -2 and T — 10 K, 
the times to reach runaway growth of the first core ranges from 0.45 Myr to 9.53 Myr 
(see Table 2). Since our simulations start from a flat density background, these times 
represent upper limits to the time that fragmentation might take for each set of param- 
eters. However, an advantage of the uniform background density is that it allows for a 
self-consistent modeling of the entire core formation process, without questions about the 
origi n of initially peaked density distributions used in earlier axisymmetric calculations 
(e.g. ICiolek fc Mouschoviasl Il993t iBasu fc Mouschoviaslll994[ ). 

In a medium with initial nonlinear perturbations, the time scales for all sets of param- 
eters are indeed likely to be shorter. However, we believe that our calculated time scales 
are relevant if the corresponding dimensional values are obtained from higher starting 
values of column density brought about in certain regions by pre-existing (including tur- 
bulent) flows. For example, the Taurus molecular cloud has an overall background number 
column density N « (1 — 2) x 10 21 cm~ 2 but also contains embedded dark clouds with 
N 
of 



5 x 10 cm , such as HC1 2 and L1495, with in which the r e are small clusters 



10 — 20 YSO's and als o many dense cores (see Gomez et al. , 19931 : Onishi et al 



20021 : 1 Goldsmith et all 120081 ) . Our periodic model may be applied to such dark clouds 



that may themselves have been brought together by nonlinear flows originating in the 
larger cloud, external triggers, or an earlier phase of gravitational fragmentation. A sec- 
ond example application of our periodic model may be wit hin the L1688 dark cloud in 
Ophiuchus, which has N « 10 22 cm" 2 (|Motte et ah . 1998). The mean spacing of frag- 
ments in these two regions varies (| Andre et al.l . l2000h . with core edges mea sured to be at 



radii < 5000 AU in L1688 but at < 20000 AU in the Taurus dark clouds (jAndre et al 
2000). While some of the difference in spacing may be due to the different background 



column densities, other important aspects of spatial and kinematic structure may also 
arise due to different values of fio and -Pext (and f n i : o to a lesser extent), as demonstrated 
in this paper. 

Unlike a uniform density three-dimensional medium, our thin sheet actually has a 
preferred scale for gravitational fragmentation with a unique value for any given set 
of initial dimensionless parameters. Since our simulation region is always safely larger 
than this fragmentation scale, it is unlikely that the size of our system (in the x- and y- 
directions) influences the final outcome, as measured by fragment spacings, time scales to 
runaway collapse, and core mass distributions, for example. This is supported by our tests 
with runs at quadruple the size of the standard simulations. Incidentally, this is not the 
case for three-dimensional periodic box simulations, in which the fastest growing mode of 
gravitational instability is always that of the box size. We believe that a stratified medium 
is also a more physical and realistic starting point than a uniform three-dimensional 
medium. This is because the formation process of molecular clouds, or magnetic fields 
and nonlinear flows within it, will tend to set up compressed regions with a characteristic 
scale similar to the half-thickness Z$ c 2 / (7rG<7 n ,o) of our adopted background reference 
sheet geometry. That scale is related to the Jeans scale and effectively determines the 
preferred fragmentation scale, as modified by magnetic field strength, ionization fraction, 
and external pressure. 

The core mass distributions that we have compiled from our simulations are narrowly- 
peaked and do not match t he bro ader distributions commonly observed in cluster- forming 
regions (e.g. Motte et al. . 19981 ). What is the missing physics that can explain the dis- 
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Fig. 13. Histogram of masses from models with fio = 0.5,0.8.1.0,1.1,2.0 and f n i o = 0.2, P C xt = 0.1. A 
weighted average of core masses is taken, so that an equal simulated area is assigned to each of models 
1 — 5. The cores are the same as the ones in the first five panels of Fig. [8] The bin width is 0.1. 



crepancy? The explanation that is closest to the spirit of our models is that real molec- 
ular clouds start their lives with an inhomogeneous distribution of physical quantities, 
including the mass-to-flux ratio. We expect that the values of /i in a molecular cloud 
may fall within the observationally-established range of [0.5,2], i.e. within a factor of 
two of the critical value in each direction. In that case, we can make a simple esti- 
mate of a core mass distribution by adding up the histograms for the five models with 
[1q = [0.5,0.8.1.0,1.1,2.0] shown in Fig. [5] Since the model with = 1-1 is run fewer 
times but with a larger box size, we sample the number of cores necessary to give equal 
area weighting with the other models. The resulting histogram is presented in Fig. 1131 
The broad variation of peaks in the individual histograms seen in Fig. [8] leads to a 
smooth power-law tail in the high mass end of the composite histogram. Fig. [13] re- 
veals a slope dlagN/dlogM w —2 in this region, only somewhat ste eper than the value 
d log N/d log M ~ —1.5 measured for example bv lMotte et al. (1998). We note that this 
general mechanism, arising from an initially inhomogeneous distribution of mass-to-flux 
ratio, is an interesting new possibility for explaining the observed broad core mass dis- 
tributions. An important point is that a relatively narrow distribution of initial mass- 
to-flux ratios (factor of a few variation) can lead to a relatively broad distribution of 
fragment masses. This mechanism remains to be explored more generally and does not 
exclud e the occurrence of o ther mechanisms like competitive accretion in a more global 
model ( Bonnell et al. . 2003h . a temporal spread of core accretion lifetimes ( Mversl . 2000t 
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2001; Gammie et al 



Basu fe Jonesl. 120041), and turbulent fragmentatio n (e.g. iPadoan et al 
200llTillev fe Pudritzl . l2007l) . 



1997; Klessen. 



The richness of the physics revealed by our models of gravitational fragmentation up 
to runaway collapse of the first core, under conditions of varying mass-to-flux ratio, ion- 
ization level, and external pressure, pave the way for more extensive models in the future. 
Adding the effects of initial cloud turbulence, implementing a technique for integrating 
past the formation of the first generation of stars, and including some form of energy 
feedback from star formation, remain to be done. The addition of new and more complex 
effects will be facilitated by the fact that the thin-sheet approximation allows efficient 
calculation of the fragmentation process while retaining a high level of realism. In this 
co ntext we point out that the recent fully three-dimensional fragmentation simulation 
of iKudoh et all (|2007l ) with magnetic fields and ambipolar diffusion bears out the main 
physical results presented by BC04. 



5. Summary 

We have carried out a large number of model simulations to study the effect of initial 
mass-to-flux ratio, neutral-ion coupling, and external pressure on dense core formation 
from gravitational fragmentation of isothermal sheet-like layers that may be embedded 
within larger molecular cloud envelopes. Our simulation box is periodic in the lateral 
(x, y) directions and typically span four nonmagnetic (Jeans) fragmentation scales in each 
of these directions. The simulations reveal a wide range of outcomes, from the unique 
transcritical fragmentation mode into massive cores, to the pressure-driven fragmentation 
into dense clusters. We emphasize the following main results of the paper: 

(i) Fragmentation Spacing. The average spacings of nonlinearly developed fragments 
are generally in excel lent agreement with t he preferred fragmentation scale of linear 



perturbation theory ( Ciolek fe Basil [20061 ). although there is definite irregularity in 



any simulation, with variation of fragment spacings. Both significantly subcritical 
and highly supercritical clouds have average fragmentation scales (A) ~ 2ttZq 1 
where Z$ is the half-thickness of the background state. The transcritical (/iq ~ 1) 
models exhibit very large (super- Jeans) average fragment spacings, although there 
is evidence for nonlinear second-stage fragmentation in some cases, which makes 
the average spacing slightly smaller than predicted by linear theory. Variation of 
the ionization fraction by a factor of two above and below the standard value does 
not have a big effect on fragment spacing. However, an external pressure dominated 
sheet undergoes dramatically smaller scale fragmentation to form a dense cluster, 
(ii) Time Evolution to Runaway. The times i run for various models to reach runaway 
collapse of the first core varies significantly for models with differing initial dimcn- 
sionless mass-to- flux ratio /iq, neutral- ion coupling parameter r^o, and dimension- 
less external pressure P cx t- Values of t lun range from 0.45 Myr to 9.53 Myr, each 
scaling as (iV nj o/10 22 cm~ 2 )~ 1 (T/10 K) 1 / 2 . The supercritical clouds evolve much 
more rapidly than the critical or subcritical clouds, with the highly supercritical 
clouds evolving sa 10 times more rapidly than a critical cloud, for the typical level 
of neutral-ion coupling. A critical cloud in turn evolves more rapidly than a sub- 
critical cloud, but the variation is a factor of order unity for plausible initial values 
of yUo; for example the [1q = 0.5 model reaches runaway collapse in a time that is 1.7 
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times longer than for the fio = 1 model. In all cases, the time to runaway collapse 
is w 10T gjm when starting with small-amplitude white noise perturbations, where 
Tg, m is the growth time of the fastest growing eigenmode mode in linear perturba- 
tion theory. The quantity r gim itself varies from w Z /c s for highly supercritical 
models to w 10Zo/c s for highly subcritical models (for a typical neutral- ion cou- 
pling level). The effect of varying f n i.o (and hence the initial ionization fraction 
Xi t o) is that t run (X T n ^ q oc X[_q approximately for ambipolar-drift-driven (critical 
or subcritical) fragmentation, so that the canonical T gjm w 10Zq/c s and our calcu- 
lated i run s=s 100Zo/c s are both possibly subject to significant variation. A pressure 
dominated cloud with P cxt = 10 has t run about 5-6 times shorter than clouds with 
small external pressure but other parameters held fixed. 

Velocities in the Nonlinear Regime. Maximum infall speeds of neutrals can become 
supersonic on core scales in the supercritical clouds, but remain subsonic for critical 
or subcritical clouds. The latter is true even if the ionization fraction is reduced 
by a factor of two. The ion speeds in the cores closely follow the neutral speeds 
but are somewhat smaller, since the gravitationally-driven motion of the neutrals 
is generally opposed in the plane of the sheet by magnetic fields. 
Core Shapes. An extensive compilation of core shapes shows that the distributions 
have a peak that is near-circular for the subcritical and critical fragmentation mod- 
els. However, the cores become more elongated in the sheet for supercritical clouds, 
with a mean axis ratio in the sheet plane as 0.5. The half-thickness remains smaller 
than cither scmiminor or semimajor axis in the sheet, so that the cores are triaxial 
and preferentially flattened along the direction of the mean magnetic field. Prelim- 
inary comparison of our results with published deprojections of the observed axis 
ratios of cores shows the best agreement with our subcritical models. 
Core Mass Distributions. An extensive compilation of core masses shows that the 
peak of the distributions are related to the preferred fragmentation mass M g)m 
of linear theory. Transcritical fragmentation yields the largest peak masses while 
the highly supercritical and decidely subcritical limits have smaller (Jeans-like) 
peaks and similar distributions. That the peak mass is strongly selected is seen in 
the very sharp drops in the distributions for greater masses. This means that the 
observed relatively broad-tailed core mass distributions cannot be explained by a 
pure local gravitational fragmentation process in a medium of uniform background 
column density and mass-to-flux ratio. Cores defined as significantly supercritical 
regions within a larger subcritical common envelope also have a very narrowly- 
peaked distribution. However, a composite mass histogram that may mimic the 
effect of an inhomogeneous assortment of initial mass-to-flux ratios, does produce 
a broad core mass distribution that resembles its observed counterparts. 
Magnetic Field Line Structure. Contraction of cores within a supercritical cloud 
yields significant curvature of the magnetic field lines and very apparent hourglass 
morphologies, since contraction proceeds primarily with field-line dragging. The 
transcritical and subcritical clouds form cores through a process driven in large 
part by neutral-ion slip, and result in lesser curvature in the magnetic field. This 
holds out the hope of using the observed curvature of possible future observations 
of hourglass magnetic fields on the core scale as a proxy for measuring the ambient 
mass-to-flux ratio. 
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